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Abstract. Long time photometric monitoring programs of gravitational lens 
systems are often carried on using modest equipment. The resolution of such 
observations is limited and some of the images may remain unresolved. It may 
be still possible to find a full set of time delays from such a blended data. We 
discuss here a particular but interesting case when we have two light curves that 
both are blends. A suitable computational algorithm is developed and tested 
to work with computer-generated model light curves. Our method combines 
both blended sequences using the hypothetical time delays between the initial 
components of the light curves as free input parameters. The combined curves 
are then compared using statistical distance estimation. It occurs that using an 
assumption of equal magnification ratios between the components of the blends, 
we can indeed recover the whole set of time delays. 

Key words: cosmology: observations - gravitational lensing - methods: sta- 
tistical 

1. INTRODUCTION 

To find the time delays caused by differences in light paths of a gravitational 
lens system, we need at least some recognizable features in the observed light 
curves. As the longest delays in some systems can be hundreds of days, we need to 
have sufficiently long measurement sets. Long-time monitorings of such systems 
are feasible with telescopes of modest size and resolution. The best example of 
this kind of photometry is the long time series obtained by Schild et al. (1997). 
The constrained resolution can be also a problem for some large scale photome- 
try programs. This motivates us to investigate possibility to recover time delays 
from the data which are not fully resolved. The general scheme of the relevant 
algorithms was developed in Hirv et al. (2007). 

In the following we will focus on the case of four original images whose unre- 
solved observations result in two blended light curves. To recover the full set of 
time delays, we will use the principles of computing the dispersion spectra intro- 
duced in Pelt et al. (1994) and Pelt et al. (1996). Hirv et al. (2007) developed 
a similar method and applied it to a three-image system where two original light 
curves were blended together but the third curve was fully resolved. 

The method of dispersion spectra was singled out as a base for our algorithms 
because of its conceptual simplicity. There are many other methods available. For 
the latest see Burud et al. (2001), Koptelova et al. (2006), Vakuhk et al. (2006), 
Cuevas-Tello et al. (2006) and references therein. Some of these new methods can 
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be generalized to handle blended data. Some critical remarks about the method 
of dispersion spectra can be found in Gil-Merino et al. (2002) and response to the 
critique in Pelt et al. (2002). 

There are two simplifying implicit assumptions in our treatment below. First, 
we suppose that the data sets contain a substantial number of observations with 
sufficient time coverage, and secondly, we totally ignore possibility of distortions 
due to the effect of microlensing. Consequently, when applying the new method 
to real observational data, certain care must be exercised. 

Our paper is organized as follows. First we introduce the basic ingredients of 
the new algorithm. Then we present test data generation methods. In the next 
part we describe the results of the numerical tests and discuss various implemen- 
tation details which can be of use for the prospective users of the new method. 
The relevant software modules may be obtained from the authors. 

2. THE METHOD 

2.1. The continuous case 

Let us have a quasar image split into four components by an intervening gravi- 
tational lens. Formally we have four functions of the quasar source variability q{t): 
fk(t) = akq{t — tk), k — I, ...,4, where are the magnification coefficients and 
tk are the flight times due to different flight paths. Our observational equipment 
can supposedly record only two images as the close pairs of /i, /2 and of /s, 
are blended together due to insufhcient resolution. Thus the corresponding signals 
gi(t) and 32(^)7 that we are going to observe, are the following functions of the 
source variability q{t): 

gi{t) = aiq{t ~ti) + a2qit ~t2), (1) 
92{t) =a3q{t-t3)+a4q{t~t4). (2) 

As the spatial separation of /i and /2 is small, we may assume, that oi « 02 and 
similarly as ~ 04 for /a and /4. The amplification ratio between gi{t) and (?2(^) is 
then a w ai/a^. Let the time delay between fi{t) and /2(i) be Aa = t2— ti, and the 
time delay between the components of the second image Ah — t^ — t^. These delays 
are typically rather short due to nearby flight paths for the component images. 
As the paths of fi{t) and /3(t) differ significantly (larger spatial separation), the 
corresponding delay Ac = tj, — ti is the longest one. Now we can rewrite the 
Eqs. dl]) and ([2]) in terms of the first subimage fi{t) and relative time delays: 

gi(i) = /i(t) + /i(t-Aa), (3) 

g2{t) = h{t-Ac) + h{t-Ac-Ah). (4) 

To keep things easier to follow wc did not multiply the Eq. @ by the amplifica- 
tion ratio a. The fact, that 171 and 32 may have different baselines and amplitudes 
is taken into account in our computational (matching) algorithm. As a schematic 
example of the initial variability, the fi{t) is shown as a single-peaked function in 
Figure 1. Shifting it by delays Aa, Ah and Ac and adding the results as in the 
Eqs. ([3]) and @ we get the double peaked blends 51 (t) and g2{t) of the source 
variability. 
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Fig. 1. Graphical explanation of the method. See text for details. 

To recover all the three independent time delays Ac, Aa and Ab hidden in the 
light curves gi{t) and g2{t), we will combine the data using three trial delays Sc, 
Sa and 6b into artificial blends A{t) and B{t): 



Ait)^g,{t-6c)+giit-Sc-6b), 



(5) 



B{t)^g2{t)+g2{t-Sa). (6) 

If it happens, that Sc — Ac, 5a — Aa and Sb = Ab, the difference of A{t) and 
B{t) vanishes to zero and this is the situation we are going to search for. The 
composition of the artificial blends A(t) and B{t), when the trial delays match the 
initial delays, is also shown in Figure 1. For clarity we plotted the components of 
the artificial blends before and after adding. Blends and components that have 
the same origin are plotted using the same line type. As we can see, artificial 
blends have the same profile, when trial delays correspond to the initial ones, and 
the difference between A{t) and B{t) vanishes. This is the idea of our method in 
terms of the continuous and noise-free light curves. Next we will see, how this 
kind of construction can be used for real noisy and sampled data. (For detailed 
discussion on adding and subtracting of sampled noisy time-series see Hirv et al. 
2007). 



2.2. The sampled case 

To build the combined sums from input data sequences (time, magnitude, 
statistical weights) tn,gn, Wn, n = 1, 2, TV and the time shifted versions of them 



4 



A. Hirv, T. Eenmde, L. J. Liivamdgi, J. Pelt 



tm, g-rm Wm^m — 1,2,..., M, we form a table of triples 

^ ,9n+ 9rn,Wn^m, [I) 

using data from the sequence gi for blend A and data from the sequence g2 for 
blend B. Values Wn,m are computed as combined weights: 

^„„,^g„.„ ^"^"^ , (8) 

where Sn,m is the downweighting function: 

1 — , if \tn — tm\ ^ ^ (9^ 

0, if \tn - im| > cr 

and a is the downweighting parameter, which depends on sampling (it can be 
prefixed or chosen using trial calculations). 

Next, by varying trial delays Sc, da and db the artificial blends A and B are 
recalculated and weighted sums of squared differences between them are found: 

= min ^n.rni^^'^n + P - B^)^W^;^ 

We may call as the statistical distance between A and B. Regression coef- 
ficients a and f3 are needed, because artificial blends may have different baselines 
and magnification. They are recalculated for every set of trial delays. The com- 
bined weights for D^: 

are formed from and which are calculated using Eq. [S] for both artificial 
blends. 

By varying trial delays Sc, 6a and Sb over pre-given grids, we are searching for 
the global minimum of statistical distance which corresponds to the recovered 
time-delay system. The weights contain parameter a to be estimated, consequently 
we need an iterative scheme to get final dispersions. First we set a = 1 in Eq. (fTTj) . 
then solve linear weighted least squares equations to get estimates for a and /3. 
After that we insert a new value of a into weights and recompute. Normally this 
process converges in 3-4 steps. A similar approach was used by Hirv et al. (2007). 

2.3. Features and difficulties of the method 

There are three additional issues, which we have to bear in mind, before starting 
actual calculations. 

• Recovering the time-delay system is a degenerate problem. The mirrored 
values of short delays Aa and Ab are also valid. For a single data set we can get 
four equally correct solutions: Ac, Aa and Ab; Ac + Ab, Aa and —Ab; Ac — Aa, 
—Aa and Ab; and Ac — Aa + Ab, —Aa and —Ab. (Interchanging gi and g2 gives 
us four additional sets of solutions, where Ac is mirrored and Aa and A6 are inter- 
changed.) All the four solutions form detectable minima in the three-dimensional 
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grid of values. For finite sequences these minima may have shghtly different 
merit function values. Our method just finds formally the deepest minimum and 
corresponding time-delay system. The recovered set of time delays may be consid- 
ered real, if it shows up as a visually noticeable minimum in the two-dimensional 
slice of statistical distance values (see low- noise part of Figure 3) . Formal signif- 
icance estimation is possible using the bootstrap-type techniques and ideas from 
Pelt et al. (1996). 

• Our method does not work if |Aa| = |A6|. Both blends are then similar, 
and we can recover only the largest delay Ac using simplest "one-dimensional" 
dispersion spectra. Having a value for the long delay it is then in principle possible 
to recover the short delay (the same for both blends) from the combined data 
using the methods described in Geiger & Schneider (1996). The combining of two 
photometric series with estimated long delay allows sometimes (if microlensing 
effect is negligible) to get a data set with twice the original sampling rate. 

The case of |Ao| = |A6| may be promptly recognized from the plot of 
values - one of the four possible solutions has a characteristic distribution along 
straight line of values (see for instance Figure 5). We may also hit an arbitrary 
solution corresponding to mirrored arbitrary short delays, which shows up as a 
normal minimum in the two-dimensional plot of values. Hence the solutions 
where \6a\ « \6b\ should be handled with care. A three-dimensional plot of 
values would be useful here. 

The tests with simulated data-sets showed, that for a reasonably good sampling 
and low noise even one day differences between Aa and Ah values can be resolved. 

• The process of calculating the values in the algorithm for two blends is 
different from its analog for a clean curve and a blend. The calculation of 
involves differences of the observed data sums. In the case of a clean image and a 
blend, we have differences of original data points and combined sums. From what 
follows that total scatter of the differences in the now method is somewhat higher 
and statistical stability is lower. Consequently, the two blend method demands 
data with higher quality. 

Our method for two blends recovers the time delays correctly also for a clean 
image and a blend. Because of different sensitivity to noise, it is sensible to use 
proper method for the nature of a given problem. For unknown nature, it is worth 
trying both algorithms for the given input data. 

2.4- Data analysis 

The procedure of analyzing real data includes the following steps. At first, 
we need to find a suitable downweighting parameter a for a given data sampling. 
Next, we should verify if the noise level of our data is under the noise value our 
method can handle. And finally, we may analyze the observed curves to recover 
the time-delay system. 

We can use an interactive simulation for estimating the suitable downweighting 
parameter for real observational data. First, we generate artificial noise- free curves 
with some pregiven time delays, using the sampling of our real data. Then starting 
from small downweighting parameter (say a = 0.5), we move on towards larger 
ones and recalculate the plot of and recover the time delays for each a. In 
general, there is an optimal a for a given sampling which recovers the time delays 
correctly and produces clearest minimum on the D'^ surface. Once we have found 
the optimum, further enlargement of downweighting parameter will not improve 
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the results. It is also possible, that for a given sampling and time-delay system, 
there is no working downweighting parameter at all. Even for a correctly estimated 
value of cr the overall success of the algorithm depends on the length of the time 
series, noise level and absolute values of actual time delays. For the best results, a 
should not be larger than half of the shortest time delay we are going to recover. 
See also Hirv et al. (2007) where the estimation of the correct cr is presented in 
some detail. 

To verify if the noise level of our data is tolerably low for the method, we may 
use the same simulated curves that were used for estimating the a parameter. 
Taking these model sequences, adding gradually increasing levels of Gaussian noise 
and recovering the foreknown time delays, we can find the maximum tolerable noise 
level for our algorithm. Then we compare the signal to noise ratio (S/N) of the 
observational data and test data which had the maximum tolerable amount of 
noise. The S / N oi observational data should be higher than for the test data. 

Once we have found that further analysis of the observational data is reason- 
able, we perform the three-dimensional search for the minimum of D^. The values 
of our trial parameters 5a, Sb and Sc which correspond to the minimum of , can 
be considered as the real recovered time-delay system, if all the presented above 
issues and complications were taken into account. 

3. GENERATING TEST DATA 

The quasar source variability (j'„,t„ is simulated using simple random walk. A 
randomly chosen value of ±1.0 is assigned cumulatively to each step in the intensity 
scale. The initial time points are generated by using random step sizes from the 
interval [0.2,1.8] days. Then the quasar signal is shifted in time by delays Aa, 
A6 and Ac and blended as in the Eqs. ^ and For blending the intensities, 
linear interpolation is used. The blend 52 is multiplied by amplification ratio 
a = 0.8 to make things more realistic (the inherently important assumption of the 
method is that both components of a given blend have nearly equal magnification 
coefficients). Both blends, gi and §2, can be resampled using generated or real 
observed sequences of time points and linear interpolation. Finally, the Gaussian 
noise is added to the model observational noise. One example of the generated 
curves is shown in Figure 2. 

4. TESTING THE METHOD 

4-.1. Simulated data 

Currently we do not have observational sequences at our disposal, that are long 
enough, sampled well and have noise level our method can work with. So, to test 
the method, we had to build artificial sequences. We generated a 4300 day long 
(2740 points) noise-free dataset with random sampling which had only daylight 
caps; Ac = 420.2, Aa = 20.2, Ab — 56.5 days and a = 0.8. This set was also 
used for estimating the optimal downweighting parameter for the given sampling. 
Different levels of Gaussian noise were added to the computed curve to check our 
method's stability against noise. This set with added Gaussian noise is shown in 
Figure 2. 
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Fig. 2. The computer-generated blends gi (lower curve) and g2 (shifted up by 130 
units). Ac = 420.2, Aa = 20.2, Afe = 56.5 days and o = 0.8. The standard deviation of 
the added Gaussian noise is 5 units. 



For our test data the op- Table 1. Recovered time delays depending 
timal downweighting parameter on the added Gaussian noise, 
was a = 1.5. To find the max- 
imally tolerable level of noise, 
we performed a three-dimensional 
search for time delays, using the 
one day step size and the follow- 
ing limits for trial delays: Sc = 
370.. .470, Sa = -30. ..70, Sb = 
6... 106. The results of noise tests 
are given in Table 1 and in Fig- 
ure 3. As we can sec, the Gaus- 
sian noise with a standard devi- 
ation of two units introduces one 
day error in the estimates of Ac 
and Aa; the Gaussian noise with 
a standard deviation of six units gives us a seven day error in the estimate of Aa; 
and a noise level of 11 units makes an error in the estimate of Aa comparable to 
its original value. The signal to noise ratio was S/N — 30 in the case of standard 
deviation of six units. Hence, for a reasonably well sampled real data we should 
keep the S/N > 30 for the method to work properly. 
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Fig. 3. Vanishing of the detectable minimum on the surface due to observational 
errors. The standard deviations of the added Gaussian noise are (clockwise from upper 
left): 0.0, 2.0, 6.0 and 11.0. Values on the color key represent the log{D^) and spacing 
of the contours. The same type of color key is used in all two-dimensional plots. 

^.2. Problem, with Schild's data 

In Hirv et al. (2007) we applied the algorithm for a clean image and a blend 
to the observational data by Schild et al. (1997). As wo got then interesting 
results, we considered applying the algorithm for two blends as well. The optimal 
downweighting parameter a = 1.5 for Schild's time series was found. Next we 
found also the noise tolerance of the method for two blends using Schild's sampling. 
Having real but bad sampling, where the points to days ratio is 0.2 and large 
gaps are included, the working noise tolerance decreased to 3.0 units (standard 
deviation). The corresponding signal to noise ratio was S/N = 50. Next we 
compared the signal to noise ratio of the Schild's data and of our simulated curves. 
Unfortunately the S/N ratio for Schild's data occurred to be lower than the ratio 
for good enough simulated curves. Because of that we cannot use assumption 
about two blends for finding time delays from Schild's data. Schild's noise level 
was tolerable for the algorithm for a blend and a clean image, but not for the 
algorithm for two blends. (For the explanation see Section 2.3.) 

4.3. The \Aa\ w |A5| case 

Having well sampled data and low noise, it is still possible to get a solution for 
very close short delays. For example, we took Ac = 420.2, Aa = 56.5, Ah — 50.1 
days having a good sampling with daylight caps only and no noise. The given 
time delays were recovered correctly. The resulting plot of values is shown in 
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Figure 4. Even a one day difference between Aa and Ab is still tolerable for the 
method, but then the minimum on the surface is not very convincing indeed. 

The singular situation, where initial Ac = 420.2, Ao = 20.2, Ab = 20.2 days 
is shown in Figure 5. We can see a characteristic distribution along straight line 
of D"^ values and no minima, but this is not always the case, as was discussed in 
Section 2.3. 

5. CONCLUSION 

A method for estimating time delays from two blended light curves was devel- 
oped and tested. As the new algorithm is more sensitive to observational errors 
than is the algorithm for a clean image and a blend, wc; wctc not able to confirm 
the interesting results for real Schild's data obtained in Hirv et al. (2007). 

Although we did not have observational data good enough for analysis, all the 
steps of recovering time delays from real data were simulated as realistically as 
possible. For planning the real observations, one should repeat similar simulations 
to establish realistic limits for observational errors and sampling. 

We believe that two new algorithms and the classical method of dispersion 
spectra form a useful toolset to analyze data which will flow out from the extensive 
photometric programs planned. If enough data and sufhcient computing power 
will be available then we can set up a new kind of searching program. First, we 
select from a general database the records, where two nearby measured images are 
variable, then we apply the delay estimation schemes. If we find that two curves 
can be described as time shifted replicas of the single source curve or blends with 
proper delay structure then the follow up spectroscopy can be called upon. 
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